Direct Wolf summation of a polarizable force field for silica 
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I. INTRODUCTION 

Silica is by far the most abundant mineral in the earth's 
crust. 1 This makes it an interesting system to study 
in simulation. Additionally, SiC>2 shows a wide range 
of crystalline structures depending on temperature and 
pressure, and it can also be solidified as a glass. Al- 
though there have been enormous advances in ab initio 
simulations of silica, 2 many effects are inaccessible due 
to length and time scale restrictions of these models. For 
large-scale atomistic simulations, a high-quality model of 
the interactions, a so-called effective potential or force 
field, is essential. 

Many attempts to parameterize the interactions in sil- 
ica have been made in the past thirty years, with various 
levels of computational intensity and accuracy. Some 
of the earlier potentials are still widely used, like for 
example the potential of van Beest, Kramer, and van 
Santen (BKS), 3 a pure pair potential with fixed charges 
and short-range corrections. However, it is believed that 
many-body effects are important for correctly describ- 
ing bond angles and bond-bending vibration frequencies 
in network- forming glasses like Si02- 4 ' 5 The potential 
model of Tangney and Scandolo 6 (TS) treats the oxygen 
atoms as polarizable. The dipole moments of these atoms 
are determined self-consistently from the local electric 
field, with short-range corrections to the polarization. 7 
A more detailed description of the TS potential is given 
in Sec. II A. 

A comparison of various silica force fields showed 8 that 
the polarizable ion model of TS yields significantly better 
results for many properties compared to the BKS poten- 
tial, while still leaving room for improvement. In a recent 
study by Paramore et al., 9 attempts to map the implicit 
many-body effects in the TS model to pure pairwise in- 
teractions did not lead to an accurate potential. This 
confirms that polarization effects are indeed necessary 
for a proper description of SiC"2 . 

In all potential models discussed above, the ions carry 
some charge qi and interact with a Coulomb potential. 



This leads to the classical Madclung problem: 10 deter- 
mining the energy of a condensed system with a pair- 
wise r~ x interaction. The convergence properties of the 
resulting sum require a special treatment, and a num- 
ber of methods to evaluate the pairwise r _1 sum have 
evolved, with the Ewald method 11 as the best-known. 
There, rapid convergence for the total Coulomb energy 
of a set of N ions with charge at positions rt that are 
part of an infinite system of point charges, 

N oo 
i— 1 J^i— 1 

(where = r 3 - — Ti and r,j = \rtj\) is assured by a 
mathematical trick. Firstly, structural periodicity of lin- 
ear size L is artificially imposed on the system, and in 
the resulting expression a decomposition of unity of the 
form 

1 = erfc(Kr) + erf(«r) (2) 
is inserted. The error function is defined as 

erf(Kr) :=-?=/ dte~ t2 . (3) 

V 71 " J 
o 

The Ewald splitting parameter k controls the distribution 
of energy contributions between the two terms. Thus, 
Eq. (1) can be written as 

N N oo 

E tot = ^yyy , m [ er fc( K |r„- + nL\) 

+ erf(/t|rjj + nL\)] , 

where the sum over periodic images n is primed to indi- 
cate that the i = j term is to be omitted for n = 0. Tak- 
ing the Fourier transform of the error-function expres- 
sion only, but not of the complementary error-function 
term, one can convert the conditionally convergent total 
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energy Eq. (1) into the sum of real-space and reciprocal- 
space contributions i?* ot and E^ 1 where each of these 
converges rapidly. The downside to the Ewald summa- 
tion method is the scaling of the computational effort 
with the number of particles in the simulation box: Even 
when the balance between real- and reciprocal-space con- 
tributions controlled by k is optimized, the computa- 
tional load increases at best as 0(iV 3 / 2 ). 12 For large-scale 
simulations with millions of atoms, this is insufficient. 
Additionally, the Ewald technique is limited to periodic 
systems. In recent years, alternative simulation tech- 
niques that show better scaling properties have been de- 
veloped, among them mesh-based methods or fast multi- 
pole methods. 13 The linear scaling, however, comes with 
considerable overhead. In contrast, Wolf et al. 14 pro- 
posed a direct summation technique with linear scaling 
(O(N)) for Coulomb interactions, that can easily be im- 
plemented in standard Molecular Dynamics (MD) codes. 
This so called Wolf summation takes into account the 
physical properties of the systems under study. 

To this end, one looks at the Fourier transform of the 
error function term of Eq. (4) 



k=tO i,j 
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where the self term (n = and i = j) is now included in 
the summation and subtracted again separately. Eq. (5) 
can be rewritten as 



lfel<fc c 



where S(k), with k = \k\, is the charge structure factor 



S(k) = 



2tt 
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5^^-exp(ifcTj) 



(7) 



The charge structure factor is the Fourier transform of 
the charge-charge autocorrelation function. 

In the systems of interest here, there are no long- 
range charge fluctuations; the charges form a cold dense 
plasma, screening each other. This means that for small 
wave vectors k, the charge structure factor is also small. 
If one now chooses a sufficiently small Ewald parame- 
ter k, the reciprocal-space contribution can be neglected 
altogether. As k is linked to the real-space cut-off r c , 
however, this might require a cut-off radius which is sub- 
stantially larger than the range of traditional short-range 
interactions like in metals. 

Concurrently, Wolf et al. also motivated a continuous 
and smooth cut-off of the remaining screened Coulomb 
potential V(rjj) = q. L qj erfc(/tr,j)r^- at a cut-off radius 
r c . The authors stated that shifting the pair potential 
so that it goes to zero smoothly at r = r c is equiva- 
lent to neutralizing the surface charge in a spherically 



truncated system. The strong fluctuations in the surface 
charge with varying r c inhibit the convergence to the true 
Madelung energy with increasing r c . The combination of 
(i) shifting the potential so that it vanishes smoothly at 
the cut-off, and (ii) damping the Coulomb potential to 
reduce the required cut-off radius, but only so weakly 
that the reciprocal-space term can still be neglected, is 
called Wolf summation. 

To evaluate the TS potential with the Wolf direct sum- 
mation technique, one first has to extend the formalism 
to the treatment of dipolar interactions. How this is done 
is shown in Sec. II. We provide an estimate of the er- 
rors made by the approximation in Sec. III. The directly 
summed TS potential was implemented in the (limited 
range) MD code IMD, 15 and various observables were 
determined and compared to the original TS implemen- 
tation with full Ewald summation (Sec. IV). Finally, we 
sum up the results in Sec V, where also an outlook is 
given. 



II. WOLF SUMMATION OF DIPOLE 
CONTRIBUTIONS 



A. Tangney-Scandolo potential model 

In the TS 6 force field, there are two contributions to 
the potential energy of a system: a pairwise potential of 
Morse-Stretch form, and the electrostatic interactions be- 
tween charges and induced dipoles on the oxygen atoms. 
The dipolc moments depend on the local electric field at 
the respective atomic sites, which in turn is determined 
by the arrangement of charges and dipoles. This implies 
that a self-consistent solution must be found. 

Tangney and Scandolo propose an iterative solution 
for the dipole moments, so that the dipole moment pf 
on atom i in iteration step n is 



(8) 



where a is the polarizability of atom i and E{r{) the 
electric field at position r^, which is calculated from the 
dipole moments (and charges) in the previous iteration 
step. The short-range dipole moment pf K is the contri- 
bution induced by short-range repulsive forces between 
anions and cations, that TS included following Rowley 
et al. 7 Starting from initial electric field strengths E°(r) 
extrapolated from the previous three time steps, Eq. (8) 
is iterated until convergence is achieved for each MD time 
step. 

The parameters of the TS potential were determined 
solely from ab initio results with the Force Matching 
method. 16 There, the potential is parameterized using 
first principles values of forces, stresses and energies in 
series of reference structures. 
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B. Smooth cut-off 

For MD with limited-range interactions, the potentials 
and their first derivatives must go to zero continuously at 
a cut-off radius r c ; otherwise, atoms crossing this thresh- 
old might get unphysical kicks. For the Morse-Stretch 
pair potential, this is generally not problematic, as it de- 
cays with Tij fast enough. In MD, following Wolf et a/., 14 
the potential UMs( r ij) is replaced by 

UMs(rij) = U M s(nj) - U M s(r c ) - (n 3 - r c )^Ms( r c), (9) 

where a prime denotes a derivative with respect to r. 

The other functions used in the TS model have a gen- 
eral r dependence of the form r~",n <E {1,2,3}. Espe- 
cially the Coulomb energy with its r~ 4 dependency can- 
not simply be cut off without a treatment as in Eq. (9), 
for otherwise the energy of the system would fluctu- 
ate strongly with r c , without convergence to the proper 
value. But even with a smooth cut-off (9), with which 
the Coulomb energy does converge, a rather large cut-off 
radius would be required to make shifting of the poten- 
tial negligible. Fortunately, the Wolf direct summation 
method 14 includes a weak exponential damping of the 
Coulomb potential by erfc(Kr). Such a damped poten- 
tial can be cut off smoothly at a much smaller radius r c 
without affecting the result. All integer powers of r -1 are 
treated in a way to conserve the differential relationship 
between the functions, i.e. the damped functions are 



r- 2 = 



dr 



r 1 -»r 1 erfc(Kr) =: /_i(r), (10) 
2ftexp(-K 2 r 2 ) 



d(r x ) d(r 1 erfc(«r)) 



dr 

- r~ 2 erfc(«r) 
: f-2(r)- 



/irr 



(11) 



This procedure is also required to conserve the energy 
during an MD simulation, as discussed in more detail in 
Sec. IIC. 

The damped potentials arc then shifted to zero and 
zero derivative at the cut-off radius, as in (9). This al- 
lows for limited-range MD simulations with a standard 
MD code. The computational effort of such a simulation 
scales linearly in the number of particles (as the num- 
ber of interactions that need to be evaluated per particle 
does not increase with the number of particles), but scales 
roughly with 0(r^). 



C. Energy conservation 

In MD simulations, the energy is conserved, if the 
forces on the particles are exactly equal to the gradient 
of the potential energy with respect to the atomic coor- 
dinates. Otherwise, the energy might oscillate or even 
drift off if not controlled by a thermostat. In standard 



MD simulations, the requirement is usually automati- 
cally fulfilled: The forces are calculated as the derivative 
of the potential, which depends directly on the atomic 
positions. In the TS model, there is also an indirect de- 
pendence, as the potential is also a function of the dipolc 
moments: 



U = U({r i },{p i ({r j })}). 



(12) 



This would in principle lead to an extra contribution to 
the derivative of the potential, 



dU dU dU d{pi\ 
d{r t } d{ri\ d{p l } d{rj}' 



(13) 



which would be practically impossible to be determined 
effectively. Luckily, if the dipole moments are iterated 
until convergence is reached, we are at an extremal value 
of the potential energy, with dll/d{pi} = 0, and so this 
part need not be evaluated. Imperfections in convergence 
may lead to a drift in the energy, however, as was already 
observed by Tangney and Scandolo. 6 

When applying the Wolf formalism to the TS poten- 
tial, another issue arises concerning the conservation of 
energy. It can most easily be explained with a simple one- 
dimensional example. Given are two oppositely charged 
point charges ±q at a mutual distance r. If the nega- 
tively charged one is polarizable with polarizability a, it 
will get a dipole moment p = aq/(kr 2 ), with k = 47re - 
This leads to a total interaction energy 



U 



-2 ■ 



2k r 



2 • 



lqp_ 
2k r 2 



q-q 



from which it follows that 



dU 
dp 



k r 2 a 



q-p 



0. 



lp 2 
' 2aJ 

dipole 



(14) 



(15) 



Here, q — q denotes the Coulomb interaction between 
charges, q— p the interactions between charge and dipole, 
and the last term is the dipole energy. When we now 
damp and cut off the interactions, we replace the r _1 , r~ 2 
functions by their damped and smoothed counterparts 
f-2(r)- If energy conservation is to be main- 
tained, the differential relation between the /_„ must 
be the same as for the r~": 



df-i(r) 
dr 



-/-2(r). 



(16) 



As a consequence, the first two derivatives of the 
smoothed damped Coulomb potential must be zero at 
re- 
in MD simulation it is computationally advantageous 
to represent pair potential functions internally as func- 
tions of r 2 , and their derivative as /' := r~ 1 df/dr. The 
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damped Coulomb potentials /_i and r 1 /_2(r) in their 
smoothly cut off version become 



/_i(r 2 ) = /_ 1 (r 2 )-/_ 1 (r c 2 ) 

- |/'-i(r 2 ) . (r 2 -r c 2 ) 



2\2 



(17) 



is then iterated until convergence is achieved. The it- 
eration starts from an extrapolation of the local electric 
field at the previous three MD time steps. To improve 
the convergence of Eq. (8), E^ d is modified after each it- 
eration step n to include a small part c from the previous 
iteration, 



and 



(19) 



lf-2(r) = lf-2(r 2 ) ,f-i(r 2 ) 



i fn 



f"-i{r z ) 



(r 2 -rl). 



(18) 



In this way, Wolf summation can be applied to dipolar 
interactions in the TS potential model. In Sec. Ill we will 
discuss why this approximation is physically justified. 



D. Implementation 

The ITAP Molecular Dynamics (IMD) package 15 is 
a flexible, highly scalable MD code for limited-range 
interactions, providing linear scaling up to thousands 
of CPUs. For finite-range interactions, the number of 
potential interaction partners of an atom is uniformly 
bounded. In order to reach linear scaling in the number 
of atoms, it is essential to find these interaction part- 
ners efficiently. IMD uses a combination of link-cells and 
neighbor lists, where the former are used to compute the 
latter in an efficient way. Since Wolf summation requires 
a relatively large cut-off radius, these neighbor list can 
get fairly big, but on today's machines this is not a prob- 
lem. Parallelization is done via a fixed geometric domain 
decomposition, where each CPU gets an equal block of 
material. For the force computation, atoms at the surface 
of a block are exchanged with the neighboring CPUs. 

All potential functions used in IMD are tabulated, even 
if some of these functions may be specified by giving the 
parameters of an analytic formula. In that case, poten- 
tial tables are constructed from the analytic formula in 
a pre-processing step. During the simulation loop, the 
functions are then evaluated by table lookup and inter- 
polation. This has proven to be the most flexible and 
efficient scheme, allowing also for very complicated po- 
tential functions. For all potential functions depending 
on the radius, care is taken that they vanish smoothly at 
the cut-off radius, along with their first derivative. 

In contrast to other interactions implemented in IMD, 
the TS potential requires a self-consistency loop within 
each time step, during which the dipole strengths of the 
oxygen atoms are determined. Before entering this loop, 
the "static" contributions -E s tat to the on-site electric 
field caused by the charges of anions and cations, and 
the short-range dipole contributions pf R are calculated 
and stored. For the "induced" part of the electric field 
.Eind, which is generated by the oxygen dipoles, Eq. (8) 



This damps the self-consistency loop and thus suppresses 
overshooting the optimal solution and subsequent oscil- 
lations. For optimal performance, a value of c = 0.2 was 
used. 

Convergence is achieved, when the root mean square 
deviation of all Cartesian dipole moment components be- 
tween two iterations is less than a user-specified tolerance 
(given in units of the dipole moment). While a larger tol- 
erance will reduce the iteration steps to convergence, it 
will also introduce a larger error in the energy conserva- 
tion, which might lead to a temperature drift in micro- 
canonical simulations. In practice, a convergence limit 
smaller than 1CP 6 Ae (with elementary charge e) will 
not lead to further improvement. With this tolerance, 
about five iterations steps are typically needed per MD 
step. 

In a parallel simulation, each CPU deals with a block of 
material. For the parallel evaluation of the energies and 
forces, at each MD step the types and positions of atoms 
near the surface of a block are first communicated to the 
neighboring CPUs. Each CPU can then perform a part 
of the energy and forces computation locally. As each 
force is computed only once, certain force and energy 
contributions have then to be communicated back to the 
home CPU of the corresponding atom, where it is added 
up. This scheme is valid for all finite range interactions. 
Since only communication between neighboring CPUs is 
necessary, the scheme is highly scalable. 

For the TS potential the procedure is very similar, ex- 
cept that now there are additional data to be commu- 
nicated. In each step of the self-consistency loop for the 
induced dipoles, the electric fields and dipole moments of 
atoms at the surface must be distributed to the neighbor- 
ing CPUs, and collected again after they have been up- 
dated. There are several additional communication steps 
for each MD step, but these arc of the same kind as for 
other short-range interactions (to neighbor CPUs only), 
and the balance between communication and computa- 
tion is not affected. For this reason, simulations with the 
TS potential will scale as well as with other short-range 
potentials. 
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III. CONVERGENCE AND ERROR 
ESTIMATION 

A. Formal Analysis 

The total interaction energy of N dipole moments Pi 
at positions is given by the expression 



0.25 



E 



1 N / 1 \ 



with Tij := Ti — Tj and r^- 



(20) 



Imposing structural 



periodicity and inserting a decomposition of unity of the 
form 



1 = erfc(Kr) + erf («r) 



(2) 



where n is again the Ewald splitting parameter, we can 
rewrite above equation as 




FIG. 1: fe-dependence of the dipole structure scalar Q(k). For 
small k, the dipole structure factor is negligible. 
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( erfc(K|rjj + nL\) + erf(«;|rjj 
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Pi- (21) 



The total energy splits into a real- and a reciprocal-space 
part: 



E 



El 



TTltot 



(22) 



Since we later intend to neglect the reciprocal-space term 
for the Wolf summation, we are interested in the contri- 
bution of E^. For its fc-behavior we have to take the 
Fourier transform of 



N 



iEE^(v®V) 



l,j 77, =0 



erf (n\rij + nL\) 
Wu + nL\ 



The prime has been omitted, since the self term (for 
n = and i — j) is now finite. Because of the three- 
dimensional periodicity the above expression can be ex- 
panded into a Fourier series: 



2nNe 2 
V 



OO 

E 

fc^O 



fc*Q(fc)fc 



cxp(-fc 2 /4K 2 



(24) 



where V is the volume of the simulation cell and Q(fe) 
the dipole structure factor 



Q(fe) 



1 



N 

iV^ 2 " E 

hi 



Pi ® P 3 e 



ih-rij 



(25) 



with the normalization factor 1 / %/ Ne 2 , where e denotes 
the elementary charge. As we can see in Eq. (24), 
the large k contributions to E%' t tend to zero rapidly, 
whereas the small k contributions are governed by the 
behavior of Q(fc), which is expected to vanish as k — >• 0. 



B. Discussion 

To legitimate the neglecting of the reciprocal-space 
term for the Wolf summation we have simulated liquid 
silica with 4896 atoms, where we get no spontaneous po- 
larization as a first result. The total dipole moment is 
p = 6.76 • 10 -29 Cm, which is insignificantly small com- 
pared to a fully polarized system and thus can be taken 
as a fluctuation. All values which arc calculated in the 
course of the simulation are time-averaged over the full 
simulation time of one picosecond. 

To analyze the k — > behavior we calculated the dipole 
structure scalar, 



Q(k) = (fc*Q(fc)fc ) s , 



(26) 



where the angular brackets indicate an average over a 
spherical shell S with width Ak centered at constant 
|fe| = k. Note that for a for a periodic system Q is 
not a continuous function, but a discrete set, consisting 
of all reciprocal space vectors. Hence the average over 
the spherical shell is necessary. Fig. 1 shows the dipole 
structure scalar in liquid silica simulations. For small 
absolute values of k, Q(k) goes to zero. 

Fig. 2 shows the fc-dependence of the reciprocal-space 
term, 



E k (k) 



2TrNe 2 exp(-fc 2 /4K 2 ) 
Q{kj — 



V 



k 2 



(27) 



for different Ewald splitting parameters k (again aver- 
aged over a spherical shell). As mentioned above, due to 
the exponential damping, large-A: contributions are neg- 
ligibly small, whereas the small-fc values are governed by 
the behavior of Q(k) as k — >• 0. 

Finally the sum in Eq. (24) is evaluated for the given 
fc-mesh with truncation sphere in the reciprocal-space. 
The difference between this approach of a spherical trun- 
cation and the full summation is very small because of 
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FIG. 2: fc-dependence of the reciprocal-space term Ek(k) for 
different Ewald splitting parameters re. The k — > behav- 
ior of Ek(k) is governed by Q(k), which results in negligible 
contributions of the small fc-values to the total energy. 
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FIG. 3: Logarithmic plot of the reciprocal-space term 
for different Ewald splitting parameters re. For sufficiently 
small re, there is no noticeable contribution to the total energy 
compared to the real-space part. 



the exponential damping in Eq. (27), as seen in the rapid 
decay of Ek(k) for increasing k in Fig. 2. In Fig. 3 the 
K-dependence of the reciprocal-space term E£ ot is illus- 
trated in a logarithmic plot. For the chosen damping of 
k = 0.1 A we get 

4^fe 0t = 3-3 /ieV, (28) 



which is small compared to the real-space part and can 
thus be neglected. 



FIG. 4: Equation of state of liquid silica for damped and 
smoothly cut off TS potential compared to experiment, 17 ab 
initio simulations and classical simulations with BKS and the 
full TS potential. 6 

IV. RESULTS 

The damped and smoothly cut off TS potential was 
used to study the same thermodynamic and structural 
properties the original authors 6 examined for the Ewald- 
summed potential. 

A. Equation of State and Bonding Properties 

We compare the equation of state of liquid silica at 
3100 K to experiments, 17 ab initio results and, of course, 
the full TS potential in Fig. 4. Pressures were obtained 
as averages along constant-volume MD runs of approx- 
imately 10 ps following 10 ps of equilibration and with 
simulation cells containing 4896 atoms. We reproduced 
the good agreement of the full TS potential with the ex- 
perimental results; both the full TS potential and our 
damped and smoothly cut off TS potential match even 
better with experiment than the ab initio results. As 
already mentioned by the original authors 6 the BKS 
model systematically underestimates the volume by ~ 
13%. The large scatter of the ab initio results can be ex- 
plained with the system size and time constraints of this 
method: especially for low pressures, the system cannot 
be equilibrated completely. 

On a microscopic level, the Si-O-Si angle distribution 
was determined from multiple MD simulation runs at 
3100 K and various pressures. The results are shown 
in Fig. 5, and are in agreement with the full TS potential 
and ab initio results. 



7 




60 70 80 90 100 110 120 130 140 150 160 170 180 
Si-O-Si angle [degrees] 



FIG. 5: Oxygen centered angle distribution in liquid silica for 
the new potential compared to simulations with BKS and the 
full TS potential as well as ab initio calculations. 
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FIG. 6: Percentage of iV-fold coordinated silicon atoms in 
liquid silica at 3100 K as a function of pressure compared to 
simulations with BKS and the full TS potential. 6 

In Fig. 6 the percentage of TV-fold coordinated silicon 
atoms in liquid silica at 3100 K as a function of pressure 
is illustrated. Our results are compared to simulations 
with the full TS potential, 6 which agree rather well with 
ab initio 18 results. 

To sum up, the equation of state and the bonding prop- 
erties of liquid silica, which the original authors 6 exam- 
ined for the Ewald-summed potential, can be reproduced 
very well by using the damped and smoothly cut off TS 
potential, while using dramatically less CPU time. Due 
to the linear scaling of computational effort in the system 
size, this advantage becomes even more pronounced the 
larger the system is. 



B. Crystal Structure Data 

We also probed the damped and smoothly cut off TS 
potential by simulating the most important low pressure 



TABLE I: Quartz 





Expt." 


New Potential 


TS" 


BKS" 


a (A) 


4.916 


4.872 


4.925 


4.941 


C (A) 


5.405 


5.359 


5.386 


5.449 


P (g/cm 3 ) 


2.646 


2.718 


2.665 


2.598 


Si-O-Si (°) 


143.7 


142.1 


144.5 


148.1 



"Reference 19. 
'Reference 6. 



TABLE II: Cristobalite 





Expt." 


New Potential 


TS" 


BKS" 


a (A) 


4.969 


5.015 


4.936 


4.920 


c(A) 


6.925 


6.999 


6.847 


6.602 


P (g/cm 3 ) 


2.334 


2.268 


2.412 


2.515 


si-o-si o 


146.4 


147.1 


144.0 


143.9 



"Reference 20. 
'Reference 6. 



crystal structures quartz, cristobalite and coesite. The 
relevant equilibrium variables density, Si-O-Si angle and 
the lattice parameters at 300 K are given in Tables I, II 
and III. The average relative deviation of the data from 
the experimental results is ss 0.9%, which is a compara- 
tively good agreement. By contrast, the BKS potential 
differs by ss 2.1% on average. Note that simulations with 
the full TS potential yield a relative deviation of the pa- 
rameters that averages at merely ss 0.7%. This decrease 
in precision might be countered by redetermining the pa- 
rameters for the smoothed and damped TS force field, as 
we suggest in Sec. V. Additionally, the ordered crystals 
might be more susceptible to spontaneous polarization 
compared to the liquid, however we could not confirm 
this in our simulations. 

It should be noted, however, that the TS potential 
was optimized to reproduce atomistic properties of liquid 
Si0 2 at 3000 K. For this reason, its application to low- 
temperature crystalline systems should be closely mon- 
itored. In the case of cristobalite we found that both 
the full TS potential and the smoothly truncated poten- 
tial energetically favor a slightly different orientational 
arrangement of the fundamental Si04 tetrahedra at low 
temperatures, with only little consequence on quantities 
given in Tab. II. 



TABLE III: Coesite 





Expt. a 


New Potential 


TS" 


BKS' 


a (A) 


7.136 


7.123 


7.165 


7.138 


6(A) 


7.174 


7.161 


7.162 


7.271 


c(A) 


12.369 


12.347 


12.377 


12.493 


n°) 


120.34 


120.34 


120.31 


120.76 


P (g/cm 3 ) 


2.921 


2.940 


2.933 


2.864 


Si-O-Si (°) 


143.6 


144.2 


144.0 


150.5 



"Reference 21. 
'Reference 6. 
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V. CONCLUSION 

In this work, we have demonstrated that the advan- 
tages of the TS polarizable force field can be captured 
and reproduced in MD simulations with a strictly finite 
interaction range. To this end, we have shown that the 
Wolf summation technique, i.e. smoothly cutting off the 
damped long range real space part of the electrostatic 
interaction, and neglecting the reciprocal space part al- 
together, is justified for the TS dipolar force field for 
silica. With a suitably large real space cut-off, the errors 
in the forces and energies are acceptable for the systems 
of interest. This can also be seen in simulation results: 
Our Wolf-summed TS potential can reproduce the ex- 
perimental and ab initio structural properties of silica 
reasonably well compared to the full TS interaction. 

By omitting the reciprocal space contribution, simula- 
tions with our potential can be performed with a stan- 
dard finite-range MD code like IMD. Thus, it can profit 
from the linear scaling of computational effort with sys- 
tem size common to this method. Similarly, the calcu- 
lations can easily and efficiently be parallelized, opening 
the door to large-scale calculations impossible with the 
standard Ewald summation technique. Moreover, once 
the reciprocal space part can be neglected, there is no 
longer any need for periodic boundary conditions. It has 
been shown that Wolf summation performs very well also 
for open or mixed boundary conditions, 14 opening up a 
wealth of new possibilities. 

As a rule of thumb, the real space cut-off radius re- 



quired for Wolf summation has been estimated as about 
five times the largest nearest neighbor distance of oppo- 
site charges in the system. 22 For silica, this amounts to a 
moderate value of about 8 A. But even with a more con- 
servative choice of 10 A, for more accurate simulations, 
for a system with 4896 atoms we obtained a speedup 
of more than two orders of magnitude compared to the 
original code of Tangney and Scandolo. Also this perfor- 
mance increase makes the new method very interesting, 
and opens up new possibilities. 

The original TS potential parameters were optimized 
for the full Ewald treatment of long-range interactions. 
Redetermining the parameters for the smoothed and 
damped TS force field with the actual cutoff used in 
simulation might improve the potential further. Addi- 
tionally, using a more flexible short-range interaction 
than the Morse-Stretch potential suggested by Tangney 
and Scandolo could lead to even better results. We plan 
to implement the TS polarizable oxide potential model 
in our Force Matching code potfit 23 to perform this 
optimisation. This implementation could then be used 
to determine polarizable oxide potential parameters also 
for other materials like alumina or magnesia. 
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